Subdiffusive master equation with space dependent anomalous exponent: 'Black 

Swan' effects 
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We derive the fractional master equation with space dependent anomalous exponent. We analyze 
the asymptotic behavior of corresponding lattice model both analytically and by Monte Carlo sim- 
ulation. We show that the subdiffusive fractional equations with constant anomalous exponent /x in 
a bounded domain [0, L] are not structurally stable with respect to the non-homogeneous variations 
of parameter fi. In particular, the Gibbs-Boltzmann distribution is no longer the stationary solution 
of the fractional Fokker-Planck equation whatever the space variation of the exponent might be. 
We analyze the random distribution of ti in space and find that in the long time limit, the proba- 
bility distribution is highly intermediate in space and the behavior is completely dominated by very 
unlikely events. 

PACS numbers: 05.40.-a 
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I. INTRODUCTION 

The last decade has seen increasingly detailed develop- 
ment of the fractional equations describing the anoma- 
lous transport in physics, biology, chemistry [1-4]. Spe- 
cial attention has been paid to slow subdiffusive trans- 
port for which mean squared displacement is sublinear 
< x 2 (t) >^ , where /x is the anomalous exponent 
/i < 1. Subdiffusion is experimentally observed for pro- 
teins and lipids on cell membranes [5], RNA molecules 
in the cells @, transport in spiny dendrites Q, etc. The 
major feature of this process is the absence of charac- 
teristic microscopic time scale. The theory of anomalous 
subdiffusion leads to fractional partial differential equa- 
tions involving memory effects. If we introduce the prob- 
ability density function p(x, t) for finding the particle in 
the interval (x, x + dx) at time t, then the subdiffusive 
transport of the particles under the influence of external 
time-independent force can be described by the fractional 
Fokker-Planck (FFP) equation 



dp 



T>\ ^L FP p 



(1) 



with 



L F pp 



d(v„(x)p) , d 2 {D tl {x)p) 



dx 



dx 2 



(2) 



The Riemann-Liouville derivative T>\ is defined as 



V\^p(x,t) 



1 d 



p (x, u) du 



T(fi)dtJ (t-u) 



(3) 



and the anomalous exponent fi < 1 is assumed to be 
constant. 

The central result of this paper is that the subdiffusive 
fractional equations with constant \x in a bounded do- 
main [0, L] are not structurally stable with respect to the 
non-homogeneous variations of parameter \x. It turns out 
that the space variations of the anomalous exponent lead 
to a drastic change in asymptotic behavior of p(x, t) for 
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Fig. 1: Non-uniform distribution of anomalous exponent n(x) 
on the interval [0, L]. 



large t. To show this high sensitivity to non-homogeneous 
perturbations, one can consider the following exponent 



/j,(x) = /j, + 6v{x) 



(4) 



with constant fi and perturbation Sv(x) (see Fig. 1). The 
asymptotic long-time behavior of the density p(x, t) with 
(U]) is quite different from that of the solution to JT]) with 
the constant value of /z. It means that the standard sub- 
diffusive equation with constant [i is not a robust model 
for subdiffusive transport in heterogeneous complex me- 
dia. 

Now let us explain our main result. The standard 
way to deal with the fractional equation like (JJ in the 
bounded domain [0, L] is a method of separation of vari- 
ables [1] . Let us consider the case of the reflecting bound- 
aries at x = and x = L when |T]) has a stationary 
solution p st {x) satisfying D M (a;)p s t = d(D fJi (x)p st )/dx. 
We can write a partial solution of ((T|) as p (x, t) — 
Pst(x)Q(x)T(t), then the time evolution is described by 
fractional relaxation equation 



(5) 



where A is the separation constant. The slow relaxation 



2 



process from the initial distribution po (x) is described by 

oo 

p(x,t) = Pst (x) (-A„i M ) Q„(x)p n, (6) 



n=0 



where po n = L po(x)Q n (x)dx and Q n (x) are the eigen- 
functions of 



LppQ —\Q. 

Here the operator L* FP is the adjoint of Lfp 



j F pk. 



8 2 Q 
dx 2 



(7) 



(8) 



II. FRACTIONAL MASTER EQUATION WITH 
SPACE DEPENDENT ANOMALOUS EXPONENT 

The question is how to take into account the non- 
uniform distribution of the anomalous exponent fi. We 
cannot simply substitute the expression like @ into ([1}. 
So we need a fractional master equation with space de- 
pendent /tt(x). Chcchkin, Gorcnflo and Sokolov were the 
first to derive the fractional diffusion equation with vary- 
ing fractional exponent (To| . They studied a compos- 
ite system with only two separate regions with different 
anomalous exponents and found interesting effects in- 
volving non-trivial average drift. A similar phenomenon 
has been analyzed in terms of two equations with differ- 
ent exponent by Korabel and Barkai [TTj . 



A. Hazard function and structured probability 
density function 

Here we present an alternative derivation which is valid 
for a general space and time dependent jump densities. 
Consider a 'space-jump' random walk model in one space 
dimension. The particle movement can be described as 
follows. It waits for a random time (residence time) T x at 



(see details in 8, 9]). The only difference between stan- 
dard Fokker-Planck equation and FFP equation is the 
rate of relaxation of p(x,t) — > p s t(x). In the anoma- 
lous case the relaxation process is very slow and it is de- 
scribed by a Mittag-Leffler function (— A n t M ) with the 
power-law decay as t — > oo. The exponential decay 
exp (— X n t) is recovered for p, = 1. 

In this paper we show that if we consider nonuniform 
perturbations of anomalous exponent as (@|, this relax- 
ation picture is completely changed. The method of sep- 
aration of variables does not work for space dependent 
/i(x). The asymptotic behavior of p(x,t) as t — > oo is 
essentially different from that given by ©• It turns out 
that in the limit t — > oo the probability density p (x, t) 
concentrates around the point x, where the perturbation 
Su(x) is located, while the stationary distribution p s t(x) ^\ x i 0) = 
is completely irrelevant (see Fig. 2 and Fig. 3). 



each point x in space before making a jump to another 
point. The index x indicates that the waiting time T x 
depends on a space coordinate x. It is convenient to 
define the hazard function [l2[ as the escape rate of a 
walker from the point x 

/ \ r Pr(r < T x < t + h | T x > r) 
7 x,t =lim . (9) 

Next step is the introduction of the structured probabil- 
ity density function £(x,i, r) that the particle position 
X(t) at time t is in the interval (x,x + dx) and its resi- 
dence time T x at point x is in the interval (r, r + dr) . The 
advantage of the structured density £ is that a random 
walk can be considered as Markovian. This is a stan- 
dard way to deal with non-Markovian processes (see 
also [ili - fl5| ). This density £(x, t, r) obeys the balance 
equation 



3£ d£ 
dt dr 



-7(x,t)£. 



(10) 



Here we consider only the case when the residence time 
of random walker at t = is equal to zero, so the initial 
condition is 



£(x,0,t) =po(x)S(t), 



(11) 



where po(x) is the density for the initial position A"(0) 
The boundary condition at t = can be written as [12j 



7(x, t)£(x — z, t, t)w(z\x — z, t)drdz, 



(12) 



where w(z\x, t) is the probability density for jumps z from 
the point x at time t (jumps are independent from the 
residence time). 

Our purpose now is to derive the fractional master 
equation for the probability density 



p(x,t) 



£(x, t, r)dr. 



(13) 



It is convenient to introduce the integral escape rate 



i(x, t) 



7(r, x)£(x,t,T)dT 



and integral arrival rate 

j(x,t) =£{x,t,0) 



(14) 



(15) 



as the density of particles with zero residence time. The 
boundary condition (|12p can be rewritten as 



3 ( x i *) 



i (x — z, t) w (z\x — z, t) dz. 



(16) 



Differentiation of (jT3|) with respect to time and substitu- 
tion of d^/dt from (fTU]) together with (jTHJ) gives 

— = / i (x — z, i) w (z\x — z, t) dz — i(x, t). (17) 
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To close this equation we need to express the escape rate 
i(x, t) in terms of p(x, t). We solve (fTU)) by the method of 
characteristics 

£(x,t,T) =e(x,t-T,0)e-io T 7(x, S )d S) T<t ( 18 ) 

Here we recognize the survival function [l2| 

*(x, r) = Pr {T x > t} = e" ^ T^' s ) ds (19) 
so the structural density £ can be rewritten as 

Z(x,t i T)=j{x,t-T)*{ X ,T), T<t. (20) 

The residence time PDF cj)(x, r) is related to j(x, r) as 



(x, r) = —d^/dr = y(x, r) exp ( — / 7(2;, s)rfs 

( 21 ) 

The balance equation for p (x, t) can be found by sub- 
stitution of (fT8|) and the initial condition £(x, 0,r) = 
Po(x)<5(t) into ([11 



p{x,t)= / j(a;,u)*(x,£-u)du+po(a:)*OM)- (22) 
Jo 

To obtain the equation for i(x,t) we substitute (|18[) and 
the initial condition into (Tl4l) 



i(x,t) = / j(x,u)cf>(x,t — u)du + po(x)(j)(x,t). (23) 
Jo 

Using the Laplace transform in (|22[) and (J23J) we elimi- 
nate j(x,t) and obtain [lCl ] 



j(x,t) = / K(x,t-T)p(x,T)dT, (24) 
Jo 

where K(x, t) is the memory kernel defined by its Laplace 
transform 



K(x,s) 



(j> (x, s) 
*(x, s) 



(25) 



B. Anomalous subdiffusion in heterogeneous 
media 



Let us consider the anomalous subdiffusive case with 
the survival probability [la ] 



tf (x, t) = E Kx 



t 



r(x) 



, < (i(x) < 1, (26) 



where [z] is the Mittag-Leffler function. The Laplace 
transforms of 'J (x, t) and 4>(x, t) are 

ii,,,,. ^ ,;,,,,) 



1 + (sr(x)) 



1 + {st{x)Y {x) ' 
(27) 



The Laplace transform of the memory kernel K (x, t) is 

s l-fj,(x) 



K{x,s) 



T ( x y(x) 

and the integral escape rate i (x, i) can be written as 



(28) 



i (x, i) 



1 



p(x,i) . 



(29) 



Substitution of this expression into (|17l) gives the frac- 
tional master equation 



dp 

at 



v 



l—fj,(x—z) 



P(x- z,t) 



-(x - z y( x ~ z ') 



w (z\x — z, t) dz 



V 



l-fi{x) 



p{x,t) . 



(30) 



where T> 



1 -/»(») 



is the Riemann-Liouville fractional deriva- 



tive with varying order. This equation can be used to 
derive the general Fokker-Planck equation [17]. If we as- 
sume that the anomalous exponent p and time parameter 
t are independent from coordinate x, this equation can 
be rewritten in terms of Caputo derivative 



,d^p 
at" 



P (x — z,t) w (z\x — z, t) dz — p(x, t). (31) 



It should be noted that the fractional equation with the 
Caputo derivative cannot be served as a model for sub- 
diffusion in heterogeneous media with varying in space 
anomalous exponent /i(x). 

Master equation (13T))) can be a starting point for de- 
riving nonlinear fractional equations. If instead of p we 
consider the mean density of particles p and assume that 
jump PDF w (z) depends on p, then one can write 



dp 
Of 



v]-^p{x-z,t) 



-(x - zy^- z ) 



w (z\p(x — z, t))dz 



V 



p(x,t) . 



(32) 



Expansion of this equation in z can give a variety of frac- 
tional non-linear PDE's. As an example, let us consider 
the case of the symmetrical kernel w(z\p) for which the 
first moment J R zw (z\p(x, t) dz = 0. Then ([32]) can be 
approximated by a non-linear fractional equation 



dp a 2 [d^v, 



■.\-li{x) 



at 



dx 2 



(33) 



with varying anomalous exponent p,(x) and nonlinear 
fractional diffusion coefficient D^(p) : 

D M = ■ ™a ( P ) = I z 2 w (z\p) dz. (34) 



First, let us consider random walk on a lattice with 
the space size a. We denote the probability of a particle 
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moving right and left from the point x as r(x) and l(x) 
correspondingly (r(x) + l(x) = 1). Then the jump pdf 
can be written as 

w (z\x) — r(x)S(z — a) + l(x)6(z + a). (35) 

The fractional master equation ([30)) takes the form 



dp 

~dt 



r(x — a) 



-V 



l—fj,(x — a) 



t(x - a)^-a) * 



p (x — a, t) + 



t(x + a y L{ - x+a ) 
i 

- n -t 



-(x)^( x ) 



p(x + a,t) — 
v]-^p(x,t). (36) 



In the limit of small a and t(x) jl8| one can obtain from 
(I36p the FFP equation with varying anomalous exponent 



dp _ d(v,(x)vl~^p) d^(D,{x)v]-^ 
dt ~ 



Ox 



dx 2 



(37) 

with the finite values of the fractional diffusion coefficient 
D^(x) and fractional drift v^ix) : 

(38) 

Note that in order to keep the fractional drift w M (a;) finite 
as a — > 0, we need to assume that r(x) — Z(x) = 0(a). 

If we put the reflecting barriers at x = and x = L 
and consider constant exponent and diffusion D^, then 
the FFP equation (|3"T)) admits the stationary solution in 
the form of the Gibbs-Boltzmann distribution 

1 f x 

p st {x) = Cexp [-U(x)\ , U(x)=-—J v fl (z)dz 

(39) 

with C^ 1 = J Q exp [— [/(a;)] dx. 

If is constant, the fractional time derivative does not 
affect the Gibbs-Boltzmann distribution 0, [2l| . But this 
result is structurally unstable with respect to any non- 
uniform variations of \x. Let us show now that the Gibbs- 
Boltzmann distribution (|39|) is absolutely irrelevant for 
the long time behavior of the solution to the FFP equa- 
tion (|37|) with non- uniform distribution of fj,(x) (j4j. 



C. Discrete model 

We divide the interval [0,L] into n discrete states. At 
each state i, the probability of jumping in the neigh- 
borhood to the left or right is given respectively by li 
and ri (li + ri = 1). The fractional equation ((36]) for 
Pi(t) = Pr{X(t) = i} can be rewritten as 



i—m-i 



Pt-i(t) U+iV 



l pi+l(t) 



Ti_i«- 



'"Pi® 



(40) 



subject to the conditions l\ — r_i = 0, r% = 1 and l n = 1, 
r n = l n +i = 0. Note that the FFP equation (|3"7) is just 
a continuous approximation of Eq. (|40[) . Taking the 
Laplace transform of (|4"0)) and using X^P*( S ) = s' we 
obtain 



SPz(s) 1 + 



i+1 



1 



(sTi. 



+ 



i+1 



(ST i+ i)^ + i 



0"*- 1 (st 4+ i)^+i (sny 

- E s ft-( 6i i : - 1411 



If one jUjvf is smaller than the others (/jm < /i,; Vi), one 
can find that spi(s) — »■ and spm(s) — »■ 1 as s — > 0. It 
means that in the limit i — > oo, we obtain 



Pi(i) -4 0, p M (i) -4 1. 



(42) 



This result in a continuous case can be rewritten as 
p (x,t) —> S(x — x min ) as t — > oo, where x mm is the point 
on the interval [0, L] at which fi(x) takes its minimum 
value. A similar result was obtained for a symmetrical 
random walk in [15j in the context of chemotaxis (anoma- 
lous aggregation). Note that Shushin [19J considered a 
two-state anomalous system with different anomalous ex- 
ponent fj, and found that in the long time limit the prob- 
ability is located in the slower state (see also [ill \2ty). 



III. MONTE CARLO SIMULATIONS 

To validate our results, we run Monte Carlo simu- 
lations with the following procedure. Random num- 
bers with uniform distribution, u and v, are gener- 
ated and then transformed into Mittag-Leffler distributed 
random numbers using the following inversion formula 

sin(/^7r) 



„ - -rlog(u) ( t s a ' n n ( ^) ~ cos(/i7r)j * [22j (see for de- 
tails We take L = 1 and divide the interval [0, 1] 
into 100 subintervals. We use r« = 1/2 + 5a(l — 2a«)/2, 
1 < i < 100 and a — 1/100. This corresponds to 



r(x) = l/2 + 5o(l/2-x), 



(43) 



so the drift v fl (x) — 10(1 — 2x)D f2 and the potential 
U(x) = 5(1 — 2x) 2 /2. All the random walkers start in 
the same state i = 40, their number N = 10 4 , n = 10~ 4 
for all i, and the long time limit is set at T = 10 5 . 

First step is to compute the exact stationary PDF 
given by (f39|) and see how well our Monte Carlo sim- 
ulations work. Fig. 2 shows that the Monte Carlo simu- 
lations agree with the Gibbs-Boltzmann distribution. 

The next step is to show that the Gibbs-Boltzmann dis- 
tribution (|39l) is absolutely irrelevant as far as the long 
time behavior of non-uniform system is concerned. The 
anomalous exponent [ii is assumed to be 0.5 for all states 
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Fig. 2: Long time limit of the solution to the system (|40| 
with fii — 0.5 for all i. Gibbs-Boltzmann distribution is rep 
resented by the line. 
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Fig. 3: Long time limit of the system (|40l 
to a perturbation. The parameters are fii 
i = 42 for which jU42 = 0.3. 



80 



when im is subject 
= 0.5 for all i except 



except one, i = 42, for which ^42 = 0.3. One can see from 
Fig. 3 that in the long time limit the probability is con- 
centrated at state i = 42. One can conclude that there 
is a complete break down in the predictions based on the 
FFP equation with uniform anomalous exponent. If the 
system was structurally stable we would expect to see 
something more like Fig. 2 again. However, the outcome 
is completely dominated by the perturbation (142 = 0.3. 
This result has a huge implication for modelling anoma- 
lous subdiffusivc transport of proteins, porous media, etc. 
In reality the environment in which anomalous transport 
takes place is never homogeneous. 

Several attempts have been made to take into account 
the random distribution of anomalous exponent (see, for 
example, [2J, |25|). One can introduce PDF f (fi) for a 
random /1 and write down the distributed-order fractional 
FPE as 



^d^p 
dP 1 



f (p) dfj, = L FP p. 



(44) 



Let us show that if we generate the random field /i(x) 
along the space interval [0,1], the asymptotic behavior 
of p{x, t) will be quite different from that of the average 
fractional equation flU]) . 

Fig. 4 shows the PDF / (/i) which will be used to 
generate the discrete uncorrelated random field /ij. The 
probability is concentrated around the point 0.6 such that 



Fig. 4: The PDF / (/i) of random anomalous exponent \x. 
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Fig. 5: One sample of the discrete random field /m along i for 
1< i < 100. 



Pr{0.5 < fi < 0.7} = 0.98. This distribution is chosen so 
that extreme values are highly unlikely to occur, with a 
purpose to show that the extreme low values dominate 
the long time behavior. Fig. 5 shows one sample of 
random field fi(x) on the interval [0, 1] which is subdi- 
vided into 100 subintervals (1 < i < 100). Fig. 5 shows 
clearly that the values of Hi fluctuate around the mean. 
The value at fj,g2 = 0.01245 has a very small probability, 
since Pr{/i < 0.02} = 2.5 x 10~ 4 . It is a very unlikely 
event, yet one can see from Fig. 6 the state i = 82 com- 
pletely dominates the long time outcome of (|4T)|) . This 
phenomenon can be interpreted as a 'Black Swan'. The 
distribution of p(x, t) is highly intermediate for large t, 
so the average behavior described by (|44l) can be very 
misleading. It has been found 24[ that the distribution 
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Fig. 6: Long time limit of the system 
random field represented in Fig. 5. 



40|) when fii is the 
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of the anomalous exponent in Eq. (|44|) leads to ultra-slow 
kinetics, but the stationary distribution is still given by 
the Gibbs-Boltzmann distribution [2l[ . Our results show 
that random space variation of the anomalous exponent 
leads to completely different behavior in the long time 
limit (see Fig. 6). It should be noted that anomalous 
diffusion is just an intermediate asymptotic. When time 
tends to infinity we expect a cross-over from anomalous 
diffusion to normal diffusion, and then we will recover 
the Gibbs-Boltzmann distribution. 

The standard tool for studying a subdiffusion is a sub- 
ordination technique [261 ] with constant anomalous expo- 
nent. It would be interesting to apply similar technique 
if possible to non-homogeneous case. It would be also in- 
teresting to take into account chemical reactions together 
with non-uniform anomalous exponent [2~7j ] . 

IV. CONCLUSIONS 

We have demonstrated that when the anomalous ex- 
ponent [i depends on the space variable x, the Gibbs- 
Boltzmann distribution is not a long time limit of the 
fractional Fokker-Planck equation. Even very small vari- 



ations of the exponent lead to a drastic change of p(x, t) 
in the limit t — > oo. We have derived the fractional mas- 
ter equation with space dependent anomalous exponent. 
We analyzed asymptotic behavior of corresponding lat- 
tice model in a finite domain with n states with different 
exponents. We have found that in this situation the prob- 
abilities pi (t) do not converge to the stationary distribu- 
tion. To illustrate our ideas, we ran Monte Carlo simula- 
tions which show a complete break down in the predic- 
tions based on the FFP equation with uniform anomalous 
exponent. Further, we have shown that the idea of tak- 
ing into account the randomness of anomalous exponent 
fi by averaging the fractional equation with respect to the 
distribution /(^t) is not applicable to a non- homogeneous 
finite domain. Monte Carlo simulations show that for ev- 
ery random realization of n(x) the PDF p(x,t) is highly 
intermediate, so the average behavior can be misleading. 
Although it is possible in theory to have a completely 
homogeneous environment, in which /i is uniform, it is 
not useful in any real application like chemotaxis (l5j 
or morphogen gradient formation [2c| because any non- 
homogeneous variation destroys the predictions based on 
this model in the long time limit. 
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